####################################################################################
###
###   Immigration, Public Housing and Support for the French National Front
###   Gloria Gennaro
###
###   Paper Table A5
###
####################################################################################


library(dplyr)
library(AER)
library(multiwayvcov)

################################################################################
# Set Up
################################################################################

# Load Working Directories
source(paste0(wd_main, '/2_code/00_working_directories.R'))

# Load data and clean
source(paste0(wd_code, '/02_data_load_clean_rdd.R'))

# Keep only years where there is an election
df = df_did[which(!is.na(df_did$fn)),]

# Keep only year>=2015, where I have demand data
df = df[!is.na(df$hlm_active_dem),]


################################################################################
# Variable definitions
################################################################################

df$imm_quant = cut(df$imm_share_99,
                   breaks=c(quantile(df$imm_share_99, probs = seq(0, 1, by = 1/3))),
                   labels=c("1","2","3"), include.lowest=T)

df$hlm_fr_rate = (1 + df$hlm_active_dem_fr)/(df$hlm_active_dem_fr + df$hlm_closed_dem_fr + 1)
df$hlm_exc_rate = (1 + df$hlm_active_dem_exc)/(df$hlm_active_dem_exc + df$hlm_closed_dem_exc + 1)


################################################################################
# Functions
################################################################################


LinearCombTest <- function (lmObject, vars, .vcov = NULL) {
  if (is.null(.vcov)) .vcov <- vcov(lmObject)
  beta <- coef(lmObject)
  sumvars <- sum(beta[vars])
  se <- sum(.vcov[vars, vars]) ^ 0.5
  tscore <- sumvars / se
  pvalue <- 2 * pt(abs(tscore), lmObject$df.residual, lower.tail = FALSE)
  matrix(c(sumvars, se, tscore, pvalue), nrow = 1L,
         dimnames = list(paste0(paste0(vars, collapse = " + "), " = 0"),
                         c("Estimate", "Std. Error", "t value", "Pr(>|t|)")))
}



estimate_by_imm = function(outvar, bdw){
  temp = df[which(abs(df$running)<as.numeric(bdw)),]
  container = data.frame()
  
  print(paste0('############ ', outvar, ' ############'))

  controls = " + scale(hlm_share_99) + scale(res_owned_share_99) + scale(imm_share_99) + scale(fn_1995)"
  formula2sls = paste0(outvar, ' ~ sruT*factor(imm_quant) + running_pre', controls, ' | dummy*factor(imm_quant) + running_pre', controls)

  print(paste0('----- 2sls ----------'))
  model = ivreg(data=temp, as.formula(formula2sls))
  vcv = cluster.vcov(model, cluster = temp$agglo_name)
  model2 = coeftest(model, vcov = vcv)
  
  container = rbind(container, c(LinearCombTest(model, c("sruT"), vcv), 'Low Immigration', bdw ))
  container = rbind(container, c(LinearCombTest(model, c("sruT","sruT:factor(imm_quant)2"), vcv), 'Medium Immigration', bdw ))
  container = rbind(container, c(LinearCombTest(model, c("sruT","sruT:factor(imm_quant)3"), vcv), 'High Immigration', bdw ))
  
  colnames(container) = c('coef', 'se', 'tval', 'pval', 'group', 'bdw')
  container[,1:4] =  sapply(container[,1:4], function(x) round(as.numeric(as.character(x)), 3))
  
  output = list()
  output[[1]] = model
  output[[2]] = model2
  output[[3]] = container
  
  return(output)
}


##################################################################################
# Estimations
##################################################################################


bdw = 900

# Foreigners
contex1 = estimate_by_imm('scale(hlm_exc_rate)', bdw = bdw)
contex2 = estimate_by_imm('scale(asinh(hlm_active_dem_exc))', bdw = bdw)
contex3 = estimate_by_imm('scale(asinh(hlm_closed_dem_exc))', bdw = bdw)

# French
contfr1 = estimate_by_imm('scale(hlm_fr_rate)', bdw = bdw)
contfr2 = estimate_by_imm('scale(asinh(hlm_active_dem_fr))', bdw = bdw)
contfr3 = estimate_by_imm('scale(asinh(hlm_closed_dem_fr))', bdw = bdw)


##############################################################################
#  Regression table
##############################################################################

cleaning = function(input){
  temp = input[,] %>% as.data.frame() %>% round(3)  # Native unmet demand rate
  temp$`Std. Error` = paste0('(', temp$`Std. Error` , ')')
  temp$`Estimate` = ifelse(temp$`Pr(>|t|)`<0.01, paste0(temp$`Estimate`, '***'),
                           ifelse(temp$`Pr(>|t|)`>=0.01 & temp$`Pr(>|t|)`<0.05, paste0(temp$`Estimate`, '**'),
                                  ifelse(temp$`Pr(>|t|)`>=0.05 & temp$`Pr(>|t|)`<0.1, paste0(temp$`Estimate`, '*'),temp$`Estimate`
                                  )))
  temp = temp[c('Estimate', 'Std. Error')]
  temp$obs = nobs(input)
  return(temp)
}

# Extract marginal effects

print('-------------- Column 1 -------------')
print(cleaning(contfr1[[2]])) # Native unmet demand rate

print('-------------- Column 2 -------------')
print(cleaning(contex1[[2]]))  # Immigrant unmet demand rate

print('-------------- Column 3 -------------')
print(cleaning(contfr2[[2]]))  # Native pending

print('-------------- Column 4 -------------')
print(cleaning(contex2[[2]])) # Immigrant unmet demand rate

print('-------------- Column 5 -------------')
print(cleaning(contfr3[[2]]))  # Native satisfied

print('-------------- Column 6 -------------')
print(cleaning(contex3[[2]]))  # Immigrant unmet demand rate
